#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"
#include <math.h>
#include <time.h>

double f(double x, double y)
{
  return (-2*(x*(x-1)+y*(y-1)));
}

double max(double x, double y)
{
  return x < y ? y : x;
}

double u_ref(double x, double y)
{
  return x*(x-1)*y*(y-1);
}

int main(int argc, char **argv) {
  if(argc < 2)
    {printf("expected: %s taille_mat epsilon nb_iteration\n",argv[0]);
      return 0;}
  else
    {
      srand(time(NULL));
      int max_ite = atoi(argv[3]);
      int N = atoi(argv[1]);
      double epsilon = strtod(argv[2],NULL);
      double h = 1/((double) N-1.0);
      // taille des blocs et nbre de blocs
      // TODO : gérer dynamiquement
      int nb_bloc;
      int size_bloc;
      
      int myrank, tag=99, tag2=98;
      MPI_Status status,status2;
      MPI_Init( NULL, NULL );
      MPI_Comm_size(MPI_COMM_WORLD,&nb_bloc);
      size_bloc = N/nb_bloc;
      
      MPI_Comm_rank( MPI_COMM_WORLD, &myrank );      
      //Création de la grille de size_block
      MPI_Comm GridComm;
      int dims[2] = {nb_bloc,1};
      int period[2] = {0,0};
      MPI_Cart_create(MPI_COMM_WORLD,2,dims,period,1,&GridComm);
      int coord[2];
      int gridrank;
      MPI_Comm_rank(GridComm,&gridrank);
      MPI_Cart_coords(GridComm,gridrank,2,coord);
 

      int i,j,k;
      double * u = malloc(sizeof(double)*N*(size_bloc+2));
      double * unew = malloc(sizeof(double)*N*(size_bloc+2));
      double * tmp = NULL;
      double maximum;
      //init u
      for (i=0 ; i<N*(size_bloc+2) ;i++)
	{
	  u[i] = 0.0;
	}
      
      int e,z;
      for (e=1 ; e<size_bloc+1 ; e++){
	for (z=0 ; z<N ; z++){
	  u[e*N+z] = (double) (rand()%10);
	}
      }
      //printf("rank = %d,u(0) = %f\n",myrank,u[N*(size_bloc+1)]);
      k=0;
      int res = 0;
      //boucle de calcul
      //printf("MPI_PROC_NULL : %d", MPI_PROC_NULL);
      do{
	//envoie de la ligne en plus
	int upperrank;
	int underrank;
	
	MPI_Cart_shift(GridComm,0,1,&upperrank,&underrank);
	
	//printf("rank = %d, upperrank = %d, underank = %d\n",gridrank,upperrank,underrank);
	
	//envoi de la premiere ligne au processus au dessus
	MPI_Sendrecv(u,N,MPI_DOUBLE,upperrank,tag,
	   &u[(size_bloc+1)*N],N,MPI_DOUBLE,underrank,tag,GridComm,&status);
	
	
	
	//envoi de la dernière ligne au processus dessous
	MPI_Sendrecv(&u[(size_bloc)*N],N,MPI_DOUBLE,underrank,tag2,
	   u,N,MPI_DOUBLE,upperrank,tag2,GridComm,&status2);
	

	
	/* if(gridrank == 1 && k==0) */
/* 	  { */
/* 	    printf("rank %d, receive %f \n",gridrank, u[0]); */
/* 	  } */
/* 	if(gridrank == 0 && k==0) */
/* 	  { */
/* 	    printf("rank %d , receive %f \n",gridrank,u[N*(size_bloc+1)]); */
/* 	  } */
	

	maximum=0.0;
	//starter
	int i_deb = coord[0]*size_bloc;
	int j_deb = 0;
	//printf("gridrank = %d, cord[0] = %d\n", gridrank,coord[0]);
	//calcul
	for (i=1 ; i<size_bloc+1 ;i++)
	  {
	    for (j=1 ; j<N-1 ;j++)
	      {
		unew[i*N+j] = 0.25*(u[(i-1)*N+j]+u[(i+1)*N+j]+u[i*N+j-1]+u[i*N+j+1]-h*h*f((i-1+i_deb)*h,(j+j_deb)*h));
		maximum=max(maximum,fabs(unew[i*N+j]-u[i*N+j]));
	      }
	  }
	//recopie
	tmp = u;
	u = unew;
	unew = tmp;
	k++;
	//transfert et mise à jour max/epsilon.
	int loc_cond = ((k>max_ite) || (maximum<=epsilon));
	MPI_Allreduce(&loc_cond,&res,1,MPI_INT,MPI_LAND,GridComm);
	MPI_Barrier(GridComm);
	if(k%1 == 0)
	  {
	    printf("rank = %d, k = %d, maximum = %f, res = %d, loc_cond = %d\n",gridrank,k,maximum,res,loc_cond);
	  }
	
	
	
      }while(!res);
      
      
      //pas de gather, on compare en paralléle.
      int m;
      int l;
      double max_diff = 0.0;
      if(gridrank == 0)
	{
	  FILE * fd = fopen("first.data","w+");
	  for(m=1 ; m<size_bloc; m++)
	    {
	      for(l=0 ; l<N; l++)
		{
		  max_diff = max(max_diff,fabs(u[m*N+l]-u_ref((l+coord[0])*h,m*h)));
		  fprintf(fd,"%f %f \n",u[m*N+l],u_ref((l+coord[0])*h,m*h));
		}
	    }
	}
      if(gridrank == 1)
	{
	  FILE * fd = fopen("Second.data","w+");
	  for(m=1 ; m<size_bloc; m++)
	    {
	      for(l=0 ; l<N; l++)
		{
		  max_diff = max(max_diff,fabs(u[m*N+l]-u_ref((l+coord[0])*h,m*h)));
		  fprintf(fd,"%f %f \n",u[m*N+l],u_ref((l+coord[0])*h,m*h));
		}
	    }
	}
      //Reduce erreur pour affichage.
      double erreur;
      MPI_Reduce(&max_diff,&erreur,1,MPI_DOUBLE,MPI_MAX,0,GridComm);
      if (myrank == 0)
	{
	  //difference solution
	  printf("Différence max avec fonction solution : %f\n",erreur);
	}
      free(u);
      free(unew);

      }
  
  MPI_Finalize();
  return 0;
}

